{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>发生时间</th>\n",
       "      <th>开关机状态</th>\n",
       "      <th>加热中</th>\n",
       "      <th>保温中</th>\n",
       "      <th>实际温度</th>\n",
       "      <th>热水量</th>\n",
       "      <th>水流量</th>\n",
       "      <th>加热剩余时间</th>\n",
       "      <th>当前设置温度</th>\n",
       "      <th>用水停顿时间间隔</th>\n",
       "      <th>事件编号</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>18742</th>\n",
       "      <td>2014-11-10 22:00:38</td>\n",
       "      <td>开</td>\n",
       "      <td>开</td>\n",
       "      <td>关</td>\n",
       "      <td>37°C</td>\n",
       "      <td>25%</td>\n",
       "      <td>26</td>\n",
       "      <td>17分钟</td>\n",
       "      <td>50°C</td>\n",
       "      <td>0.033333</td>\n",
       "      <td>170</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18743</th>\n",
       "      <td>2014-11-10 22:00:42</td>\n",
       "      <td>开</td>\n",
       "      <td>开</td>\n",
       "      <td>关</td>\n",
       "      <td>37°C</td>\n",
       "      <td>25%</td>\n",
       "      <td>23</td>\n",
       "      <td>17分钟</td>\n",
       "      <td>50°C</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>170</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18744</th>\n",
       "      <td>2014-11-10 22:00:46</td>\n",
       "      <td>开</td>\n",
       "      <td>开</td>\n",
       "      <td>关</td>\n",
       "      <td>37°C</td>\n",
       "      <td>25%</td>\n",
       "      <td>25</td>\n",
       "      <td>17分钟</td>\n",
       "      <td>50°C</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>170</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18798</th>\n",
       "      <td>2014-11-10 22:19:43</td>\n",
       "      <td>开</td>\n",
       "      <td>关</td>\n",
       "      <td>开</td>\n",
       "      <td>50°C</td>\n",
       "      <td>50%</td>\n",
       "      <td>8</td>\n",
       "      <td>0分钟</td>\n",
       "      <td>50°C</td>\n",
       "      <td>18.950000</td>\n",
       "      <td>171</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18800</th>\n",
       "      <td>2014-11-10 22:49:07</td>\n",
       "      <td>开</td>\n",
       "      <td>关</td>\n",
       "      <td>开</td>\n",
       "      <td>50°C</td>\n",
       "      <td>50%</td>\n",
       "      <td>8</td>\n",
       "      <td>0分钟</td>\n",
       "      <td>50°C</td>\n",
       "      <td>29.400000</td>\n",
       "      <td>172</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     发生时间 开关机状态 加热中 保温中  实际温度  热水量  水流量 加热剩余时间 当前设置温度  \\\n",
       "18742 2014-11-10 22:00:38     开   开   关  37°C  25%   26   17分钟   50°C   \n",
       "18743 2014-11-10 22:00:42     开   开   关  37°C  25%   23   17分钟   50°C   \n",
       "18744 2014-11-10 22:00:46     开   开   关  37°C  25%   25   17分钟   50°C   \n",
       "18798 2014-11-10 22:19:43     开   关   开  50°C  50%    8    0分钟   50°C   \n",
       "18800 2014-11-10 22:49:07     开   关   开  50°C  50%    8    0分钟   50°C   \n",
       "\n",
       "        用水停顿时间间隔  事件编号  \n",
       "18742   0.033333   170  \n",
       "18743   0.066667   170  \n",
       "18744   0.066667   170  \n",
       "18798  18.950000   171  \n",
       "18800  29.400000   172  "
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# dataExchange_thresholdOptimization.py\n",
    "# -*- utf-8 -*-\n",
    "# 用水事件阈值寻优模型\n",
    "# 第一步：确定阈值的变化与划分得到的事件个数关系\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from pandas import DataFrame\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "inputfile = 'dataExchange_divideEvent.xlsx'\n",
    "data = pd.read_excel(inputfile)\n",
    "data.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data.drop([u'事件编号'],axis=1,inplace = True)\n",
    "data.to_excel('thresholdOptimization.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[197,\n",
       " 194,\n",
       " 191,\n",
       " 186,\n",
       " 181,\n",
       " 178,\n",
       " 174,\n",
       " 172,\n",
       " 172,\n",
       " 171,\n",
       " 171,\n",
       " 171,\n",
       " 169,\n",
       " 164,\n",
       " 164,\n",
       " 163,\n",
       " 161,\n",
       " 159,\n",
       " 158,\n",
       " 158,\n",
       " 158,\n",
       " 155,\n",
       " 154,\n",
       " 153]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#确定阈值与事件数的关系\n",
    "#****************************\n",
    "#@1  目标：确定阈值的变化与划分得到的事件个数关系\n",
    "#    方法：通过频率分布直方图\n",
    "#****************************\n",
    "timedeltalist = np.arange(2.25,8.25,0.25)\n",
    "# 从2.25到8.25间，以间隔为0.25，确定阈值即，阈值范围为[2.25,2.5,2.75,3,...,7.75,8]\n",
    "counts = [] # 记录不同阈值下的事件个数\n",
    "for i in range(len(timedeltalist)):\n",
    "    threshold = pd.Timedelta(minutes = timedeltalist[i])#阈值为四分钟\n",
    "    d = data[u'发生时间'].diff() > threshold #  # 相邻时间做差分，比较是否大于阈值\n",
    "    data[u'事件编号'] = d.cumsum() + 1 # 通过累积求和的方式为事件编号\n",
    "    temp = data[u'事件编号'].max()\n",
    "    counts.append(temp)\n",
    "\n",
    "counts  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAF0CAYAAAA+UXBRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8lOW5//HPRQJJIOyUnQrIooBNjamIDYsrWrV1X6Bq\nT/UoaOvS1qq1VYv7Uo9ad6vn+HOrC9TdugNGA21AUUGRCqgoiyKURQgQrt8f95MmxkAWMvPMTL7v\n1yuvzDPzzDPX3BO45t7N3REREZHM0SLuAERERKRpKbmLiIhkGCV3ERGRDKPkLiIikmGU3EVERDKM\nkruIiEiGUXIXERHJMEruIiIiGUbJXUREJMNkxx1AQ3Tp0sX79u0bdxgiIiJJMWvWrC/d/TsNfV5a\nJfe+fftSVlYWdxgiIiJJYWYfN+Z5apYXERHJMEruIiIiGUbJXUREJMMouYuIiGQYJXcREZEMo+Qu\nIiKSYZTcRUREMoySu4iISIZRchcREckwSu4iIiIZJr2S+/z5sGxZ3FGIiIiktPRK7uvWwaRJcUch\nIiKS0tIruQPcfjuYQV5e3JGIiIikpPRL7i1bwvjxsGhR3JGIiIikpPRL7ps3w5Il0L173JGIiIik\npPRK7rvsAt/9LkybBg88EHc0IiIiKSm9knubNvDBB7DPPnDyyTB5ctwRiYiIpJz0Su4QBtI99RTs\ntReccAI8+2zcEYmIiKSU9EvuAPn58Nxz8L3vwVFHwSuvxB2RiIhIykjP5A7Qvj288AIUFIRBdiIi\nIgJAdtwB7JDOnaG0FFpE31FWr4YOHeKNSUREJGbpW3OvVJnYH3oIBgwIA+5ERESasfRP7pVGjICf\n/AR22inuSERERGKVOcm9Xz+4554wmn71aq1gJyIizVbmJPfqjjsuzIX/5JO4IxEREUm6zEzuV10V\nau/77QdLl8YdjYiISFJlZnIvLITnnw+Jff/94Ysv4o5IREQkaTIzuUMYYPfMM7BwIRx4IKxaFXdE\nIiIiSZG5yR1gzBj4299g7lw4+GBYuzbuiERERBIus5M7wEEHwaOPQlkZHHoofP113BGJiIgkVOYn\nd4DDDw9bxL7+Olx3XeiLHz0ali2LOzIREZEml97LzzbE8cdDly4wahSccw6UlMCkSXDbbXFHJiIi\n0qTM3eOOod6Kioq8rKys8RfIy4ONG799f24ubNjQ+OuKiIgkgJnNcveihj6veTTLV1q4EMaNg1at\nwnFuLowfr9XsREQkozSv5N6jB7RrB1u2QE4ObNoUjrt3jzsyERGRJtO8kjvA8uUwYQLMnFn1+7zz\nII26J0RERLan+QyoqzRlStXtW2+Fs8+G668P/fGTJsUXl4iISBNpfsm9pv/5H1i/Hi67DFq3hgsu\niDsiERGRHaLk3qIF3HlnGC1/4YUhwZ91VtxRiYiINJqSO0BWFtx3X0jwZ58dmuj/+7/jjkpERKRR\nmt+Aum3JzoaHHw5r0J9+eljRTkREJA0puVeXkwOTJ4cNZ372s3BbREQkzSi515SXB089BfvtB506\nxR2NiIhIg6nPvTb5+fD3v4NZOF62TAvdiIhI2lDNfVsqE/uDD8KAAfD22/HGIyIiUk9K7nXZf/8w\ncn7XXeOOREREpF6U3OvSrVtY6CYnB1auhLlz445IRERku5TcG+LEE2GffeCDD+KOREREZJuU3Bvi\nxhvDinb77QcffRR3NCIiIrVScm+IQYPg5ZehvDwk+E8+iTsiERGRb1Fyb6hhw+DFF2H16pDg334b\nRo8O0+VERERSgJJ7YxQWwvPPw9KlsO++UFKi7WJFRCRlKLk31r77hq1iV62CrVvh9tvD3Pi8vLgj\nExGRZk7JvbEWLoRx48IUOYCWLWH8eFi0KN64RESk2VNyb6wePaBdO9i8OST2iopwrGVqRUQkZkru\nO2L5cpgwAf75z/D788/jjkhEREQbx+yQKVOqbp98Mvz4xzBjBuy1V3wxiYhIs6eae1PZdVf44Q+h\nTZu4IxERkWYuYcndzLqZ2evR7X5m9qyZvW5mf6p2zj1m9qaZ/T5RcSRN27YweTLstlvckYiISDOX\nkORuZh2B+4DKauw1wGXuPhLobWZjzOxIIMvd9wZ6mtnARMSSdMuWwTnnwLp1cUciIiLNVKL63CuA\n44Ano+NBwOzo9gqgPTAGeDS671WgGFhQ80JmdhpwGkDPnj2ZOnUqAP3796dt27bMmTMHgM6dOzN0\n6FCmT58OQHZ2NsXFxcyePZs1a9YAUFRUxPLly/n0008BGDhwIDk5Obz33nsAdO3alUGDBlFSUgJA\nTk4OI0aMoKysjHVRsh4+fDhLlizhs88+A2Dw4MFkZWUxb948APouXUrfm27i45UrWXTKKeTl5TF8\n+HBmzpzJhg0bABgxYgSLFi1iWbSq3ZAhQ6ioqGD+/PkA9OrVi969ezNz5kwA8vPzKSoqorS0lPLy\ncgCKi4v58MMPWbFiBQDDhg2jvLycBQtCEfbp04du3bpRVlYGQLt27SgsLKSkpIQtW7YAMGrUKObO\nncvKlSsBKCgoYO3atSxcuDC8l7596dSpE7Nnh4+uY8eOFBQUMG3aNNwdM2P06NHMmTOHVatWAVBY\nWMhXX33F4sWLU/pz6t69O/369aO0tBRAn5M+J31O+pxS8nNqLHP3Rj+5zoubTXX3MVGzex4wA7gB\n2B24CbjZ3eeY2YFAobtfvb3rFRUVeeWHltLGjw9N9B98AH37xh2NiIikKTOb5e5FDX1eUgbUufvl\nwPPAqcB97r4OWEdI+AD5yYolKa65BrKy4Lzz4o5ERESaoWQm1LeB7xJq7gCzCE3xAAXA4iTGkli9\ne8P558Pjj8O0aXFHIyIizUwyk/t5wA3u/nV0/ARwopndABwLPJvEWBLvN7+B7343DK6rqIg7GhER\naUYSmtzdfUy125e4+/3VjtcQBtXNAPZx938nMpaka90arr02bAl7771xRyMiIs1IrP3c7r7K3R91\n98zcDP3YY6G4GC66CP6dWd9dREQkdWXOILZUZAY33girV8Orr8YdjYiINBNaWz7R9tgDPv447CIn\nIiKSBKq5J0NlYp87N944RESkWVByT5annoJhw+Cll+KOREREMpySe7KMHQvXXQd77x13JCIikuHU\n554sOTlh7ruIiEiCqeaebG+8AcOHQ7RpgYiISFNTck+29u2hrAwuuSTuSEREJEMpuSfbsGEwYQLc\ncQdEWy6KiIg0JSX3OEyaBO3awbnnQgK33BURkeZJyT0OnTvDH/8IL78MTz8ddzQiIpJhlNzjMmEC\n7Lor/OpXUF4edzQiIpJBlNzj0rIl/M//wEcfwc03xx2NiIhkECX3OI0dC4ceCpddBsuXxx2NiIhk\nCCX3uP3pT9CzJ3zySdyRiIhIhtAKdXEbNAjmzYMW+p4lIiJNQxklFbRoARs2wD33aGqciIjsMCX3\nVPHYY3DqqfDEEzB6NCxbFndEIiKSppTcU8X48fDmm2FL2JKSsNCNiIhII5inUTNwUVGRl5WVxR1G\nYuTlwcaN374/Nzc02YuISLNjZrPcvaihz1PNPVUsXAjjxkGrVuE4NzfU5hctijcuERFJO0ruqaJH\nj7De/ObN4XjjxnDcvXu8cYmISNpRck8ly5fDxIlhYB3A3LnxxiMiImlJ89xTyZQp4fe6dfDcc/D1\n17B1q+bAi4hIgyhrpKL8fLj6aigrg/vvjzsaERFJM0ruqWr8eBg+HC64ANaujTsaERFJI0ruqapF\nC7jpprCYzVVXxR2NiIikESX3VDZ8OJx4YthcZuHCuKMREZE0oQF1qe6qqyA7OyxyIyIiUg9K7qmu\nVy+49964oxARkTSiZvl08e67cPrpUFERdyQiIpLilNzTxQcfwOTJMH9+3JGIiEiKU3JPF0cfDR99\nBEOGxB2JiIikOCX3dGEG7duHZvnS0rijERGRFKbknm6uuAJGjQrN9CIiIrVQck83EyZA69bw61/H\nHYmIiKQoJfd007UrXHxx2Fjm+efjjkZERFKQkns6+uUvYeBAOPfcqv3fRUREIkru6ahVK7jhhjAt\n7tZb445GRERSjJJ7ujrkEBg7Fv74R/jyy7ijERGRFKLknq7MQu197drQBy8iIhJRck9nQ4bAmWfC\nnXfC++/HHY2IiKQIbRyT7i69FHbfHQYPjjsSERFJEUru6a5jR/jZz8Jt99BcLyIizZqa5TPFX/8K\ne+4JGzfGHYmIiMRMyT1TdOkCHTrAqlVxRyIiIjFTs3ym2H9/2G8/NcuLiIhq7hnFDD75JIyeFxGR\nZkvJPdPcfnvYXKasLO5IREQkJkrumebCC6FbNzj77DB6XkREmh0l90zTrl3Y8/3NN0MtfvRoWLYs\n7qhERCSJlNwz0c9+BoWF8NvfQkkJTJoUd0QiIpJESu6ZKD8fZs+G9eth69ZQgzeDvLy4IxMRkSRQ\ncs9ECxfCuHGQlRWOc3Jg/HhYtCjeuEREJCmU3DNRjx6h771yOdryclizBrp3jzsyERFJAiX3TLV8\neZgSN3VqWL1u6dK4IxIRkSTRCnWZasqUqttLl0J29FF//TW0bh1PTCIikhSquTcHlYn9+eehf394\n77144xERkYRKWHI3s25m9np0u6OZPWdmr5vZHdXOucfM3jSz3ycqDqlm112huBh69Yo7EhERSaCE\nJHcz6wjcB7SJ7joReMDdRwJtzazIzI4Estx9b6CnmQ1MRCxSTd++8PjjYQ/4jRvDOvQiIpJxEtXn\nXgEcBzwZHa8EBptZB6AP8AlwEvBo9PirQDGwoOaFzOw04DSAnj17MnXqVAD69+9P27ZtmTNnDgCd\nO3dm6NChTJ8+HYDs7GyKi4uZPXs2a9asAaCoqIjly5fz6aefAjBw4EBycnJ4L2qm7tq1K4MGDaKk\npASAnJwcRowYQVlZGevWrQNg+PDhLFmyhM8++wyAwYMHk5WVxbx58wDo3r07/fr1o7S0FIC8vDyG\nDx/OzJkz2bBhAwAjRoxg0aJFLItWjhsyZAgVFRXMnz8fgF69etG7d29mzpwJQH5+PkVFRZSWllJe\nXg5AcXExH374IStWrABg2LBhlJeXs2BBKMI+ffrQrVs3yqI15tu1a0dhYSElJSVs2bKFIZMm8Z2P\nPmL+XXexLCcHgIKCAtauXcvChQsB6Nu3L506dWL27NkAdOzYkYKCAqZNm4a7Y2aMHj2aOXPmsCra\narawsJCvvvqKxYsX63Nqgs8JYNSoUcydO5eVK1fqc9LnpM+pmX1OjWWewPXHzWyqu48xs52Aq4AP\ngN7AmcAdwM3uPsfMDgQK3f3q7V2vqKjIy7QhStN46y3YZ58wkn76dOjZM+6IRESkBjOb5e5FDX1e\nsgbUXQlMcPdJhAT/X8A6oHLJtPwkxiIAu+8Of/97mDK3//7wxRdxRyQiIk0kWQm1NbCbmWUBwwEH\nZhGa4gEKgMVJikUq7bUXPPNMWLnugAMgamISEZH0lqzkfhVwF/BvoBPwMPAEcKKZ3QAcCzybpFik\nutGj4Ykn4P334aCDwkp2IiKS1hKa3N19TPT7H+4+1N3z3f0Ad1/n7muAMcAMYB93/3ciY5HtGDsW\nHnssbDZz6KFhoRsREUlbsfZzu/sqd3/U3bXheNx+/GN44AF44w24+OK4oxERkR2g5WelynHHhe1i\nR42KOxIREdkBGqEu33TIIdC2LaxbB1deCdHcUBERSR9K7lK7p58OzfMzZsQdiYiINJCa5aV2J5wQ\n5sLvskvckYiISAOp5i7bVpnYn34aLrgAPv88TJ1bpvGPIiKpTMld6vbaa3DNNWFEfUkJTJoUd0Qi\nIrIdda4tb2b5wK+AnQkbwgA84u4vJDi2b9Ha8jHJywu7yNWUmwvRphAiItL0ErK2vJkNAP4P+Ju7\nn+zuPyds+lJkZtc3KlJJPwsXhj74FtGfS24ujB8flq0VEZGUU1ez/M+Bn7v7u5V3uPsGd78CmGdm\n30todJIaevSA9u2hspVn40Zo0wa6d483LhERqdV2k7u7/y5aJra2x+5193cSE5aknOXLYeJEuOmm\ncPzKK/HGIyIi21TvqXBmNhr42N0XJy4cSVlTplTdXrAAbrklbBl70EHxxSQiIrVqyGj5a4AhiQpE\n0si118KwYfDf/w2bNsUdjYiI1FCvmruZnQSscPfnEhyPpIO8PHjkkdD33qpV3NGIiEgNdSZ3MxtF\nGCE/NvHhSNoYUq0RZ9Ei6NcvvlhEROQbtpnczSwbuAoYDJwOHGFmFdVOaQG0cve7EhuipLQ774Sz\nzgp7wQ8dGnc0IiLC9mvubYDewDrC4jXtgOormRjQMnGhSVo46ihYsQIGDow7EhERiWwzubv7v4ET\nzOxI4Dbgx+6+KmmRSXro0gX+8Idw++uvoXXreOMREZG6R8u7+xTgeuARM7PEhyRpacECGDz4m1Pm\nREQkFvWaCufuTwKfAhMTG46krZ12CivWnXoqLFkSdzQiIs1aQ+a5XwL8PVGBSJpr1QoefjjMez/x\nRKioqPs5IiKSEPVO7u6+xN0XJjIYSXMDBsCtt8LUqWGLWBERiUVdu8JdZ2YdtvHYGdo4Rr7lpJPg\n+OPh4oth5sy4oxERaZbqqrn/GfiLme1VeYeZtTWzK4Ee2jhGvsUMbr8devcO28SuqXXfIRERSaC6\ndoX7BDgRKDaz/2dm9xBGzv/d3f+QjAAlDXXoAA89BB9/DGeeGXc0IiLNTp3Lz7r7BkJCF6m/vfeG\nSy4JP+ecA3vsEXdEIiLNRkNGy4s0zO9+B9OnK7GLiCSZkrskTnY2jBwZbs+eDZs3xxuPiEgzoeQu\niffBB7DnnnDDDXFHIiLSLNRrP3eRHbLLLnDHHXDMMXFHIiLSLCi5S3Kcemr4XV4OGzdC+/bxxiMi\nksHqbJa3oGgbj/Vr+pAkY1VUhD74n/8c3OOORkQkY9W3z/1JM7vBzH5hZj8AMLM9gIcSF5pknKws\nOPbYsHPc9dfD6NGwbFncUYmIZJy6lp81d3fgX8CdwIfAwWY2C7gFUCeqNMyvfgX77w8XXgivvw6T\nJsUdkYhIxqmr5v68mT0OtAcGAHsBRcDTwBLgO4kNTzJOmzbw8suhid49LFVrBnl5cUcmIpIx6kru\nRwOXAqXAJGBP4Ah3vxQ4D/iTmVkiA5QMs3AhjBsHublV9w0eDPPmxReTiEiGqSu5nwGcBawG3gEu\nAhaZ2VWETWXOj5rtReqnRw9o1y7s+16Z4OfPh0MPhbfeijc2EZEMUVdybxf93hvIIUydmwO8DfQF\n5iYsMslcy5fDhAkwYwaccUZYh371anj66bgjExHJCHXNc3+JkNiPB14BjgAKgC+Ba4HfAb9PZICS\ngaZMqbp9663h91dfhRo9QElJqOHvvHPyYxMRyQB11dz3BTYRRsYPBu4njJyf4e73AzubmZawlR3X\nqVNYi37rVjjtNDjxRM2FFxFppO3W3N39EjNrQxgtXxGd/3t3fzM6ZaK7b01wjNKctGgBL7wA69eH\nUfRr1sCGDdCtW9yRiYikjTpr3e6+3t0/d/flwDJgaLXHVicyOGmm+vQJ69EDnHceDBsGTzwRb0wi\nImmkXk3qZvYagLtXEKbHiSTH2WeHZH/EEWHZ2jVr4o5IRCTl1be/vHrnZ0UiAhGp1ZAhYVT9RRfB\nffdBQUFY2U5ERLapruVnS8zsJaDAzF6MbhdV3jazGckJU5q1Vq3g8stDUm/RIqxJf/75YYc5ERH5\nlu0md3cvdvcDgHfc/cDo9qzK2+6+V3LCFCHMh58zJ2wfe+21sOee8O67sHRp025C09TXExFJsrpq\n7gPNbBrQ2cxaJykmkW3Lz4e77goL3ixbBrfdBpddFubGN9UmNE19PRGRJLPtrR5rZr0JS88eDkwE\nxgH/5+77JCe8byoqKvKysrI4XlpSUW5u7U3zubmhGf+f/9z+89u1C18UAP70p9DUX1HLkJLc3DAd\nT0QkycxslrsXNfR525znbmZ5wP8CnwAfE1arOxl408wujk5rCfzR3bc0PGSRHbRoEfzmN2Ga3Ndf\nh/74E04Ie8Vffjm8/fb2n9+p0zevtd9+0KVL1fUgHJ9zTlgyV3PtRSRNbDO5u/sG4AAzKwSuBtYB\nN1A1ct4IyV2j5yUelZvQbNwYatebNoXj7t3hllsadq3K8ydODNfLyQnX27oVfv97uPhi2Hff8OXh\nyCOhQ4emfz8iIk2krj73rsA6dz8QmAYsAtYAX7p7ibu/pl3hJFbVN6GZMGHHB8FVXm/mzJDoR4+G\n996DCy8M29WeckqowR9+OMya1TTvQUSkidXV574vsBvwFLAZyCP0u68HnnT3BckIspL63CVW7qEf\n/69/hUcegcmTYa+9YO7c0Kw/diy0bBl3lCKSQRrb517XIjYVwBbgdeBu4FZgODASuLmhLyaS1szC\n9LsbboBPPoHhw8P9d98Nxx0XmvEBliypGpiXytP0NOVPJGPVd4W694EZhH72vwL/D8hKVFAiKS8r\nKyR7CHPuS0qgTZtwfNhhYcncc86BX/widafpacqfSMaqq1l+EtAd2BP4ANgFmA6UAwcAx7r7h0mI\nE1CzvKQBd3j8cTj++DAYrzZduoS++6uvhs2boWdP+N3v4Nxz4dNPobDw28/58svar5WTA23bwk03\nwbhxodvgRz/afozbupam/ImknCafChfJjX6ygZzody6hBp9NqMFrlTqRSmZwzDFQXAxnnQVPPRWa\n67OzoV8/+OEPoXVr2H33cH6LFnDssbDrruG4detwXNP69fDmm6Fvf8uWkNSPPjp8Kbj11nBtgM6d\na39+zWt9/jm88UbVlL8jjggLAolIRqgruT8LDCPU3t8C2gAlhAF1uwDXJzQ6kXTVo0eooW/ZUjVN\nb//9v51As7JCcq7UufM3j6ubODEsulN92t+QId88v3//bT+/5rVeeaVqyl/37uFn06awlr+IpLX6\n9rkPBfaueae7v9C04YhkkERN02uK69Wc8rdsGbz1FgwYAM8+u2Nxikjs6jMVbhjwKLAVGE+ota8F\nerr7n5IRZCX1uYsk0Ny5od/+nXfgtNPCkrz5+XFHJdKsJWoq3FzgeXdf5u4r3P1/gL8D05Od2EUk\nwYYOhX/8A3772zC97/vfh9LSuKMSkUaoa8vX5TUXqnH3T9z9s7oubGbdzOz16PZEM5sa/bxtZndG\n999jZm+a2e935E2ISBPJyYFrroGpU8N4geLisPxu5Rx+EUkL9e1zbxAz6wjcRxiAh7vf7u5j3H0M\n0YI4ZnYkkOXuewM9zWxgImIRkUYYNSo0z590ElxxBYwYAfPmxR2ViNRTXaPlG6sCOA54svqdZtYL\n6O7uZWZ2EqEvH+BVoBj41nK2ZnYacBpAz549mTp1KgD9+/enbdu2zJkzB4DOnTszdOhQpk+fDkB2\ndjbFxcXMnj2bNWvWAFBUVMTy5cv59NNPARg4cCA5OTm89957AHTt2pVBgwZRUlICQE5ODiNGjKCs\nrIx169YBMHz4cJYsWcJnn4XGi8GDB5OVlcW86D++7t27069fP0qj5sy8vDyGDx/OzJkz2RDNIR4x\nYgSLFi1iWTQgasiQIVRUVDB//nwAevXqRe/evZk5cyYA+fn5FBUVUVpaSnm0xWlxcTEffvghK1as\nAGDYsGGUl5ezYEEowj59+tCtWzcqxyi0a9eOwsJCSkpK2LIlbOI3atQo5s6dy8qVKwEoKChg7dq1\nLFy4EIC+ffvSqVMnZs+eDUDHjh0pKChg2rRpuDtmxujRo5kzZw6rVq0CoLCwkK+++orFixfrc8qE\nz+nkk+nSvz+D//QnWu63H29NmcK/o/emzymFPif9e8rYz6mxtjugbkeZ2dSotl55fCXwsru/amb3\nADe7+xwzOxAodPert3c9DagTicmyZaHmvu++YXGeL77QFrgiSZCoAXVNxsxaAPsCr0V3rSNsRAOQ\nn8xYRKSBuncPiR3g9tthl11gQVL3jRKRBkhmQh0JzKi2RewsQlM8QAGwOImxiEhjjR0Lp58e5sRD\nqMlrExqRlJKoPvfajCWsS1/pCeB1M+sJHIyWsRVJDwMGhHXxAT7+OKxl379/1SY0WsZWJHYJ7XOv\n88XDqPoDCPPm6/zKrz53kRRTuXxtTdqERqRJpHyfe23cfZW7P1qfxC4iKWjx4rBRTXbUCNiyJYwf\nHza4EZHYaBCbiDRejx7QqVPod2/RImxh26pVGIAnIrFRcheRHVO5Cc2j0bIV06bFG4+IJHVAnYhk\noilTqm6fcQbceWfYhGbo0PhiEmnmVHMXkabzxz9C27Zw7rkQ42BdkeZOyV1Emk6XLiHBv/EGfPRR\n3NGINFtK7iLStCZOhA8/rFrkRkSSTsldRJpWy5bQq1dollftXSQWSu4ikhgXXwy77w7RDl4ikjwa\nLS8iiXHiidCzZ5gHLyJJpeQuIokxaFD4EZGkU7O8iCTWQw/B0UdrapxIEim5i0hirVsHkyfDY4/F\nHYlIs6HkLiKJdcopUFAA552nneJEkkTJXUQSKysLbrwRPvkErr8+7mhEmgUldxFJvDFj4Kij4Oqr\nYcmSuKMRyXhK7iKSHNddBxUVcMEFcUcikvGU3EUkOfr1g1//Gh58EEpL445GJKMpuYtI8lx4IfTo\nAWefDVu3xh2NSMZScheR5MnPh2uuCYn9iy/ijkYkYym5i0hyjR8P//gHdOsWdyQiGUvJXUSSq0WL\n8LNyJTz9dNzRiGQkJXcRicfvfgcnnACrV8cdiUjGUXIXkXhccgnMnAkdOsQdiUjG0a5wIhKPnj3D\nD8D69dCmTbzxiGQQ1dxFJF6/+Q388IdhgRsRaRJK7iISr+HDYc4cuOeeuCMRyRhK7iISr6OPhpEj\n4aKLNLhOpIkouYtIvMzgppvC1LjLLos7GpGMoOQuIvHbffew7/vNN8P8+XFHI5L2lNxFJDVcfjm0\nbh02lxGRHaLkLiKpoVs3+MMf4Nln4aGHYPRoWLZsx6+7dGnTXUskTSi5i0jqOOssGDgQzjwTSkpg\n0qQdv+ZllzXdtUTShLl73DHUW1FRkZeVlcUdhogkSl4ebNz47ftzcsL9TzwRBt/V5e67Ybfdar9W\nbi5s2LANF9pzAAAX2klEQVTjsYokgZnNcveihj5PNXcRSR0LF8K4caHvHUJS79oVXnstHLuH7WLr\n+qm81siRYZMaCNccMgT23x8mT6498YtkCC0/KyKpo0cPaNcuJN7cXNi0CY46CkaMCI8fcUT4qa+h\nQ+GNN8K1Nm4MCf6f/4RnnoG2bcO1TjgB9tsPWrZMzHsSiYFq7iKSWpYvhwkTYMaM8HtHBsLVvFaf\nPrBkCbz4Ylg858kn4eCDwxr3Z5wBr79eVfMXSWPqcxeR5qu8HJ5/Hv76V3jqqdAXf9FFYVpe5f+N\nZvHGKM2a+txFRBoqJwcOPzwk9xUr4MEHQ58/wLRpMHgwzJ37zedoap2kAfW5i4gA5OdXJXaA7GwY\nMAD69g3H998fEvu771ZNrbvttlhCFamLmuVFROojO7v2bWk1tU4SSM3yIiKJ9Omn8OMfhyQPoS/+\nwANh0aJ44xKphZK7iEh99OgRRtVv3QqtWoUBdy++CLfcEqbsiaQQJXcRkfqqnFr3j3+EXey++124\n4oowD//99+OOTuQ/1OcuIrIj/vY3OO00GDasaiU9kSbS2D53jZYXEdkRRxwRau6Vy9muXBkG2PXu\nHW9c0qypWV5EZEd17141Ze6cc6CoCNavjzUkad5UcxcRaUqXXgo/+hG0aROOv/66aiMckSRRzV1E\npCntvHPYjAZCf/zAgWFUvUgSKbmLiCTKTjtB+/Ywdiz88pehFi+SBEruIiKJUlgIs2bB2WeH+fB7\n7AGa8SNJoOQuIpJIeXlw443w0kuwdm0YWX/ZZbBlizahkYRRchcRSYb99w+bzhxzDFx8MRQXw69/\nXbUJjUgT0iI2IiLJ1qoVbN787fu1CY3UoI1jRETSxccfh33kc3PDcV4eHHmkNqGRJqPkLiKSbD16\nhIVvNm2qqq0/+yy0bBl3ZJIhlNxFROJQuQnNjBlhXvygQdC5c3hsy5Z4Y5O0pxXqRETiMGVK1e2H\nHqq6XVoKP/0p3HtvGEkv0giquYuIpJLsbGjRAvbZB847D8rL445I0lDCkruZdTOz12vcd5uZHVbt\n+B4ze9PMfp+oOERE0soPfgBvvw2nnw7XXx+O33kn7qgkzSQkuZtZR+A+oE21+0YC3d396ej4SCDL\n3fcGeprZwETEIiKSdtq0gdtvD4PsVqwIu8xdey1UVMQdmaSJRPW5VwDHAU8CmFlL4G7gOTP7ibs/\nCYwBHo3OfxUoBhbUvJCZnQacBtCzZ0+mTp0KQP/+/Wnbti1z5swBoHPnzgwdOpTp06eHN5adTXFx\nMbNnz2bNmjUAFBUVsXz5cj799FMABg4cSE5ODu+99x4AXbt2ZdCgQZSUlACQk5PDiBEjKCsrY926\ndQAMHz6cJUuW8NlnnwEwePBgsrKymDdvHgDdu3enX79+lJaWApCXl8fw4cOZOXMmG6L5qyNGjGDR\nokUsi1alGjJkCBUVFcyfPx+AXr160bt3b2bOnAlAfn4+RUVFlJaWUh410RUXF/Phhx+yYsUKAIYN\nG0Z5eTkLFoQi7NOnD926daNyXYB27dpRWFhISUkJW6LBOqNGjWLu3LmsXLkSgIKCAtauXcvChQsB\n6Nu3L506dWL27NkAdOzYkYKCAqZNm4a7Y2aMHj2aOXPmsGrVKgAKCwv56quvWLx4sT4nfU76nHb0\ncyouZvVzz5F3zjl85/zz2Th5MuV33cVb0eP6nFLkc0rgv6fGSugiNmY21d3HmNkpwCHAGcAvgWXA\n94Gb3X2OmR0IFLr71du7nhaxEZFmyR3uvx9+8Qv4zndg/vzQNy8ZL9UXsdkduMvdlwEPAPsA64C8\n6PH8JMYiIpJezOCkk8LytffcExJ7RQVEtU+RmpKVUP8F9I9uFwEfA7MITfEABcDiJMUiIpKedtoJ\nxowJt2+4AYYMgSVLYg1JUlOy2nXuAe41s+OBlsDRwFrgdTPrCRwM7JWkWERE0t/YsaHm3qtXOP78\n87AYziOPhNXvpFmLdeOYaFT9AcD0qMl+u9TnLiJSiwULYK+9YNWqsOrdbbfFHZE0kcb2uWtXOBGR\ndJaXBxs3fvv+nJza75e0kuoD6kREJBEWLoRx40KShzD4DsL69IccAg88AGvXxhefxELJXUQknfXo\nAe3ahWVqc3NDcj/mGPj1r+G99+DEE6FrVzjjjLgjlSRSchcRSXfVd5ibMCHU2q+5JuwPX1ICp5wS\n5scDbN0K55wD6uLMaOpzFxFpTv71L9hjD7j11rD73PLlYUDe3nuHDWsqLV0Kxx/fNKPvm/Jaibhe\nClOfu4iI1G3AgJDQjzkmHD/wAIwcCX37hl3oZs8OK+Jddlmo9U+atOOv2ZTXSsT1MpBq7iIizdm6\ndfDkk/Dww/DCC6FJvza5ubBhA7zxRlgdb9SocP8rr2x7wN6xx8Lmzd++v1Wrqq1sn3229nOq69kT\n9txz2zMDKmPLQJoKJyIiO2blSrj33rDVbLQ5CxB2qfvXv0IT+AEHwNdfhyQPMHQoRBvIbFOLFqGv\nv3XrsHTuqFHw9NPhsY4dYfXq7T//6KPhscdCc/xOO4WWhS1bwvWOOCLEm6HN841N7tp5QEREgs6d\nQ9P8woVw553QsmWoVR96aFXyvP32kFwrTZ68/fn0V1wRzsnNDecdeyxceWXV49On172Vbfv24XeP\nHnD44fD441W19datMzax7wgldxER+ably2HiRDjtNLjrrlBjrjRgwDfP3WWX7V+rouLb1+rXr+rx\n3XZrWGxbtoTrHXssHHhg6FL485/Doj3yH2qWFxGR9HTPPXDqqaE2/+ijoaUhw2i0vIiINC+nnAI3\n3wxPPAEnn1x3834zomZ5ERFJX7/8ZRjgd8EFoR/+L3/55nz9ZkrJXURE0tv554cEP2lSGGD35z9X\nrbHfTCm5i4hI+rv00pDgr78+JPhrrmnWCV7JXURE0p8ZXHttSPDXXQdHHhn2uG+m1DEhIiKZwSw0\nyb/2WrNO7KDkLiIimaRFCxgzJtyeNg3uuCPWcOKi5C4iIpnprrtCTX57K+hlKPW5i4hIZrr33rCp\nTW5u3JEknWruIiKSmXJyoEsX2LQJxo2Dv/0t7oiSRsldREQy26ZNsGgRHHccPP983NEkhZK7iIhk\ntvz8kNSHDQtT5F57Le6IEk7JXUREMl+HDvDii7DzznDYYfDmm3FHlFBK7iIi0jx06QIvvww9e8LB\nB8OsWXFHlDBK7iIi0nx07w6vvAIdO4b94F99FUaPhmXL4o6sSSm5i4hI89KnT0jqublwyCHw+uth\n05kMouQuIiLNz9Ch8PnnYYEbd7j99rB8bV5e3JE1CSV3ERFpfhYuDHPfW7cOx3l5YV78j35UdU5F\nRTyxNQEldxERaX569IB27ULNPTcXysvhgAPgwgvD4wsWhHPOPBNKSmDr1njjbSAldxERaZ6WL4cJ\nE2DGjPC7ZUsoKgqPbdkC++wD//u/MHIk9O0L550Hs2eHZvwUZ54GQVYqKirysrKyuMMQEZHmYt06\nePJJePhheOGFkPQHD4bjj4cTTgi3q1u6NDz2yCNhZP4OMrNZ7l7U0Oep5i4iIrIt+fkwfjw880yY\nLnfnnaG5ftKkUMuv3HFuw4bw+7LLQjN+zKPvVXMXERFpqM8+g3feCYvhQNhHvrZ8mptblfgbQTV3\nERGRZOnVqyqxb9kCv/0tFBdXjb5v0QL22gvmzIklPCV3ERGRHZGdDVdfHTam2bgRWrUKo+tnzIDv\nfS9sVvPYYztUg28oJXcREZGmUDn6/h//gIkTw7K2EyZAaSkceyx07Qo//Sk8+yxs3pzQUNTnLiIi\nkkgVFTBtWhhxP3kyrFoVEv2iRVXN+NugPncREZFUlJUF++4Ld98dRtw/9RScfXZVYj/pJLjuum8/\nb+lSdoXB336gbtk7Eq+IiIg0QKtWYT/5ww4LxxUV8PXXVVPqNm2Cq66Co4+GW2+lNeQ35mXULC8i\nIpIqSkth773/c1gElLlbQy+jZnkREZFUMWIEvP027LFHaM5vJCV3ERGRVFJQAD/4Abjj0KjmdSV3\nERGRVBNNq5sP7zfm6RpQJyIikmqmTAFg/W23NWrlG9XcRUREMoySu4iISIZRchcREckwSu4iIiIZ\nRsldREQkwyi5i4iIZBgldxERkQyj5C4iIpJhlNxFREQyjJK7iIhIhlFyFxERyTBptZ+7mX0BfNxE\nl+sCfNlE12pqiq3hUjUuUGyNpdgaLlXjAsXWWIPdvW1Dn5RWG8e4+3ea6lpmVubuRU11vaak2Bou\nVeMCxdZYiq3hUjUuUGyNZWZljXmemuVFREQyjJK7iIhIhmnOyf2uuAPYDsXWcKkaFyi2xlJsDZeq\ncYFia6xGxZZWA+pERESkbs255i4iIpKRlNxjZGadzOwAM+sSdyzVpWpcIiJSPxmd3M2svZk9b2Yv\nmdnfzKxVLedkm9knZjY1+tktSbH1AJ4F9gReM7Nap/mZ2T1m9qaZ/T5V4oqrzKq9fjcze2s7jye1\nzGq89jZji/FvrV6vG0e51Se2FPh7u83MDtvO47H8vW0vrhj/1iZWe823zezObZwXx99anbHFWG4d\nzew5M3vdzO7Yznn1LreMTu7AeOAGdz8AWAYcVMs53wMedvcx0c+7SYptKHCuu18BvAAU1jzBzI4E\nstx9b6CnmQ1MhbiIr8wqXQ/k1fZATGVW3TZjI75yq/N1Yyy3+pRJbH9vZjYS6O7uT2/j8VjKra64\niKnM3P32ytcEXgfurnlOXGVWn9iI72/tROABdx8JtDWzb825b2i5ZXRyd/fb3P2l6PA7wIpaTtsL\nOMLMSszsQTNLysI+7v6yu88ws1GEWnJpLaeNAR6Nbr8KFKdIXLGUGYCZ7QusJ3xZq80YklxmleoR\nW1zlVp/XHUM85Vaf2GIpNzNrSUgAi83sJ9s4bQxJLrd6xhXbv1EAM+tF+PJR2wIsY4jp3yjUGVtc\n5bYSGGxmHYA+wCe1nDOGBpRbRif3SmY2Aujo7jNqefifwGh3LwZWAz9KYlwGHAdsBipqOaUN8Fl0\new3QLUXiiqXMLHSrXAxcsJ3T4iqz+sQW199afV43lnKrZ2xxldtJwDzgWmBPM/tlLefEUW71iSu2\n/9ciZwK3b+OxuP7WKm0vtrjKrQQYCJwFfACsquWcBpVbxid3M+sE/Bn4+TZOecfdl0a3PyAUcFJ4\ncCbwJnBoLaeso6qJN58kfV71iCuuMrsAuNXdV2/nnFjKjPrFFle51ed14yq3+sQWV7ntDtzl7suA\nB4B9ajknjnKrT1yx/b9mZi2AfYHXtnFKXH9r9YktrnK7Epjg7pOi1/2vWs5pULlldHKPalOPAhe6\n+7Y2nLnfzArMLAs4ApiTpNjON7OTosMOhG+JNc2iqumlAFicInHFUmbA/sCZZjYV+L6Z/aWWc5Je\nZpH6xBZXudXndeMqt/rEFle5/QvoH90uovZNq+Iot/rEFVeZAYwEZvi2F1GJ628N6o4trnJrDewW\nve5woLb4GlZu7p6xP8BEQvPG1OjnEuDyGucMA94B3gWuSGJsHYGXgOnAbYSBbDVja0f447oBeB9o\nnyJxxVJmNWKYCgxJhTJrQGxx/a1943VTqdzqGVtc5dYWeCz6t1AK/DAVyq2eccX2b5RQCz0yup0y\nf2v1jC2uv7U9gbmE2vlLhIF9O1RuWqEuxZlZR+AAYLqHZjipg8qscVRujaNyaziVWeM0pNyU3EVE\nRDJMRve5i4iINEdK7iIiIhlGyV1ERCTDKLmLZAAza12Pczqb2QnR7ZbRYkV1PSen+nlmlhVN18HM\nisxsLwuy6rhObrXbLaNV1rZ3ft+6YmsoM+tf91kimSGpSxKKZDoz+znwG+BzoBewFVhKmGJY4u5n\n1zj/XGCDu99R4/7uhGmIlSsEHgJsBF6JjrOAudVGzD4VvXauu39oZm+4+w9rhLcOuNrM5gIXAt3M\nbGv02B5AP+Aw4A9R3PcRVircUj2/A5cCT0fn3AucAFxrZpupqjAMAI5z98pNdJ4wsysJy2oeCIw0\nsz9E7+Njd99S7b2fDyyg2jxeM/sZ8Dhhnq+7+ws1ymsgcKO7H1LtvseAi9z9w+iuw8zsC3d/CJEM\np+Qu0rS2EDYr+ouZTQA2uvv/mdkYQoKuaSPRv8NoHeu+7v4vwpzWYVQl9+5AeXQfQEvCUpTLzGwP\n4FPCmg7PAqOADdVfJKpZdwJ+BSxz9xNqPD41un45YR7tQGAK4QvBUdVOfcWjzUrcfbaZ/Ze7zwHG\nmtk04AlgPmHe7pbo2jtH180FjgF+AOQAR0fv/VZgbXRuX6CPu19TLbadgF8Svmy8BTxjZtPcfWO1\nuC4AOpnZnwnLALcBWkX39Y7e803ReuFPu/vamh+ESCZRchdpWg6cZ2Y/BXoCW6NaZwfg72Y2GriK\nUOutTMbfj557I/AF8Meo9r2WsCrVPKAHIdH3JSy+MdLdP4qedzXwBnA4IZlNBQqi3+buownLVT7u\n7j80sxdqNKNX7pboUVwAe7j7fDM7nbAO91zCWtYnw38237gY6GBmlWtc70ZYHnMdYYGVynm2VxIW\n3XiZ0GLQO3qd9sBlNRLticAtNcr0FuB3HubtLjezB4AHzex4d99sZuMIC7vsQ/jScA3hi8XOwPGE\nDTeOAj4iLNd6OHA/IhlMyV2kaTlw3TZq7gcR1uvfl1Azr9y84idm1g542d3vqnatrYSVskoJXwA2\nERJ9BWFTH8zsHMJOdAacAYxw97Vm9rK771/tWuWVzwGy3X2/6PlT3b16szuEpvUt0YIZTxBWyxoT\nPfZgtevlEFYx/A7hS8lvCF9cWhCW48XMjiEslbnI3beaWRtCAofQktGxRvnt7O4fVB5E3RargBf/\nU8ChBt4LeMPMTiV8cehLSN4PERL7LUAXwhemDlS1gMwgdDsouUtGU3IXaVrb/Tfl7puBzVHin0xI\n8vcARTUSO8B4Qp90MVU1907RY5UbSEwm1O73A8bVbG6OBsPVHOy2i5m9HN0uqCXMrcB1wO8JX0i+\nqLwckG9mOYQvKQC/AG4iJPY1wDmEmnqludF9Y6rFPSC63bWW164e+06EL0IfAu+ZWRdgUbX3cwuh\n+X02oQn+J4Ty+Ol2LruBqrITyVhK7iJNqw0w0cyOJxpQFzXRV67Zj5l1JQxUGwmcStgFapOZnV/Z\n1xwltvWEZmQICX4jUBY9d5yZTXH3OdG5rYD/NbNyQuvB983s74RE+Dihv7rS+5W1+qjpvqYcQm13\nIKHlYAHwMHB59LtVtXOzo/u/JOwxfa67v21mhwO4+7waI/l7RO8ZwjiCl2q89gYzy3f3dR42ezos\nivMNYJi7T4re7w3ufl/0WHtgdFTGo6Pr3EvYXOWB6Pel0f39COMTRDKakrtI0xpC2LqxpLZmeQtb\nEE8G/lCjOfwS4P/M7CFCIvqaMMhsU/T4xuhnHaF5fX21xwA2ufuYyoOoWf6gase51c79fvWaezSQ\nr7p8Qh/1dYS9r38EHEn40rASqByp3pJQux9KGF/QhfAl5b8Ig+dq8yWhqR/CoLuanote+74a9/+U\n0JdO9Fr/SdDu/m/CaP0DCaP4vwf8kbDV8wWEXbbWRacfSxgoKJLRlNxFmkjUBL4XcF7lXdEPhH9r\nFcBOwN3uXplccwhT4Rw4OepDXu/un5lZT+BgQlKtbJb/fvT7SXefX+116tKi8jx371JL7K0IzfHZ\nwEp3P9bMekTn7xGNhD8AODeKJ5vQSvErwvS8T4DTovd+IiH5Vg6oM6BFNIjv30BJdH+H6LWz3L2y\nT/wZ4BEze97dV0SPnwVscfeZ0TmDqFH7jprs9wAucPeTovuI4n8kOh4E9HL3d+tRXiJpTcldpOns\nA/zT3b+OjtcQ+tf7AtcDZ0fzvt8CMLMjCU3sJ1dewN3/Uu329cD1ZnYQMAk42d3fr+V1W0Y/1bWr\ncZxP+CLxLWb2IGGQ3aYai8tY9B4ALiPUgE8HxhIS+2p3PyW6RgvCNLl/RzXofamqaecQavXPEfrv\nL632Gj8g/D/01+g9u5n9CtjfzB4m1LJXEvWjm9klhC8ZJ9V4G3mEKXyXV7svu8Z73pcwNkAk42lX\nOJEmZGat3H1T3Wc26WtmAVmNfV0za9uU875r1MR39Fpt3H19U1xLpDlRchcREckwWlteREQkwyi5\ni4iIZBgldxERkQyj5C4iIpJh/j+4AKjIGA6jqQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xa112940>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 画频率分布直方图\n",
    "#将阈值与对应的事件数绘制成频率分布直方图，以确定最优阈值\n",
    "from pandas import Series\n",
    "coun = Series(counts, index=timedeltalist)\n",
    "plt.rcParams['font.sans-serif']=['SimHei']\n",
    "plt.rcParams['axes.unicode_minus']= False\n",
    "plt.rc('figure', figsize=(8,6))\n",
    "np.set_printoptions(precision=4)\n",
    "fig = plt.figure()\n",
    "fig.set(alpha=0.2)#设置图标透明度\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "# coun.plot(linestyle='-.',color='r',marker='<')\n",
    "coun.plot(style='-.r*')#同上\n",
    "ax.locator_params('x',nbins = int(len(coun)/2)+1)  # (****)\n",
    "ax.set_xlabel(u'用水事件间隔阈值(分钟)')\n",
    "ax.set_ylabel(u'事件数（个）')\n",
    "ax.grid(axis='y',linestyle='--') # (****)\n",
    "plt.savefig('threshold_numofCase.jpg')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# 由上图可知，图像趋势平缓说明用户的停顿习惯趋于稳定，所以取该段时间开始作为阈值，既不会将短的用水时间合并，也不会将长的用水时间拆开\n",
    "# 因此，最后选取一次用水时间间隔阈值为4分钟\n",
    "# 利用阈值的斜率指标来作为某点的斜率指标"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 第二步：阈值优化\n",
    "#****************************\n",
    "#@2  目标：确定阈值的变化与划分得到的事件个数关系\n",
    "#    方法：通过图像中斜率指标\n",
    "#****************************\n",
    "\n",
    "# 当存在阈值的斜率指标 k<KS :\n",
    "#     取阈值最小的点A（可能存在多个阈值的斜率指标小于1）的横坐标x作为用水事件划分的阈值（该值是经过实验数据验证的专家阈值）\n",
    "# 当不存在阈值的斜率指标 k<KS：\n",
    "#     找所有阈值中“斜率指标最小”的阈值t1：\n",
    "#     若：该阈值t1对应的斜率指标小于KS2：\n",
    "#         则取该阈值作为用水事件划分的阈值\n",
    "#     若：该阈值t1对应的斜率指标不小于KS2\n",
    "#         则阈值取默认值——4分钟\n",
    "# 备注：\n",
    "# KS是评价斜率指标用的专家阈值1\n",
    "# KS是评价斜率指标用的专家阈值2\n",
    "\n",
    "data = pd.read_excel('thresholdOptimization.xlsx')\n",
    "n = 4 #使用以后四个点的平均斜率\n",
    "KS = 1 # 专家阈值1\n",
    "KS2 = 5 # 专家阈值2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Anaconda2\\lib\\site-packages\\ipykernel_launcher.py:9: FutureWarning: pd.rolling_mean is deprecated for Series and will be removed in a future version, replace with \n",
      "\tSeries.rolling(window=4,center=False).mean()\n",
      "  if __name__ == '__main__':\n",
      "D:\\Anaconda2\\lib\\site-packages\\ipykernel_launcher.py:12: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  if sys.path[0] == '':\n"
     ]
    }
   ],
   "source": [
    "def event_num(ts):\n",
    "    d = data[u'发生时间'].diff() > ts # 相邻时间做差分，比较是否大于阈值（*****）\n",
    "    return d.sum()+1 # 直接返回事件数（*****）\n",
    "\n",
    "dt = [pd.Timedelta(minutes = i) for i in np.arange(1,9,0.25)]#（***）\n",
    "h = DataFrame(dt,columns = [u'阈值']) # 定义阈值列（**）\n",
    "h[u'事件数'] = h[u'阈值'].apply(event_num) # 计算每个阈值对应的事件数（*****）\n",
    "h[u'斜率'] = h[u'事件数'].diff()/0.25 # 计算每个相邻点对应的斜率（****）\n",
    "h[u'斜率指标偏移前'] = pd.rolling_mean(h[u'斜率'].abs(), n)# 采用当前指标和后n个指标斜率的绝对值的平均作为当前指标的斜率\n",
    "\n",
    "h[u'斜率指标'] = np.nan\n",
    "h[u'斜率指标'][:-4] = h[u'斜率指标偏移前'] [4:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "当前时间最优时间间隔为4.0分钟\n"
     ]
    }
   ],
   "source": [
    "mink = h[u'斜率指标'][h[u'斜率指标'] < KS]# 斜率指标小于1的值的集合\n",
    "mink1 = h[u'斜率指标'][h[u'斜率指标'] < KS2]# 斜率指标小于5的值的集合\n",
    "\n",
    "if list(mink): # 斜率指标值小于1不为空时，即，存在斜率指标值小于1时\n",
    "    minky = [h[u'阈值'][i] for i in mink.index]# 取“阈值最小”的点A所对应的间隔时间作为ts\n",
    "    ts = min(minky) #取最小时间为ts\n",
    "elif list(mink1):# 当不存在斜率指标值小于1时,找所有阈值中“斜率指标最小”的阈值\n",
    "    t1 = h[u'阈值'][h[u'斜率指标'].idxmin()] #“斜率指标最小”的阈值t1\n",
    "    # ts = h[u'阈值'][h[u'斜率指标偏移前'].idxmin() - n] #等价于前一行作用（*****）\n",
    "    # 备注：用idxmin返回最小值的Index，由于rolling_mean自动计算的是前n个斜率的绝对值的平均，所以结果要平移-n，得到偏移后的各个值的斜率指标，注意：最后四个值没有斜率指标因为找不出在它以后的四个更长的值\n",
    "    if h[u'斜率指标'].min()<5:\n",
    "        ts = t1#当该阈值的斜率指标小于5，则取该阈值作为用水事件划分的阈值\n",
    "    else:\n",
    "        ts = pd.Timedelta(minutes = 4)# 当该阈值的斜率指标不小于5，则阈值取默认值——4分钟\n",
    "\n",
    "tm = ts/np.timedelta64(1, 'm')\n",
    "\n",
    "print \"当前时间最优时间间隔为%s分钟\" % tm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "172"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coun[4.00] # 可知，当最优时间是4分钟时，事件数是172个"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
